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ABSTRACT 



The complex geometry of many parts of a ceramic gas tur- 
bine engine and the need to know the stresses in these parts 
to a high degree of accuracy indicate that the finite element 
method should be employed. One analysis code used a finite 
element program based on an isoparametric 12-node brick. 

This allowed quadratic variation in one coordinate direction, 
but only linear variation in the other two coordinate direc- 
tions. This introduced an artificial stiffness into the 
structure. This study was made to determine the effects of 
using a totally quadratic isoparametric brick, vice one 
which is quadratic in only one direction. 

A distributed load preprocessor was also developed and 
tested. 

FORTRAN IV was used throughout. 
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I. INTRODUCTION 



A. WHY CERAMICS? 

There are three major advantages in using ceramics in 
the hot section of a gas turbine: improved engine perform- 

ance, higher durability, and less dependence on expensive 
imported materials . 

Simple Brayton cycle analysis shows that a higher tur- 
bine inlet temperature (TIT) results in a higher thermal 
efficiency, and lower specific fuel consumption.. Blade 
material considerations usually limit TIT. Intricate blade 
cooling schemes, such as those used in the LM2500 and the 
FT9, use impingement and convective heat transfer and exter- 
nal film cooling, to lower the temperature of the turbine 
blade. This degrades engine performance and introduces 
manufacturing complexities. Additionally, the cooling 
holes are susceptible to plugging, which may lead to even- 
tual failure. 

Ceramics, however, can withstand higher turbine inlet 
temperatures without blade cooling. Ceramic components 
also have significantly higher resistance to hot corrosion 
than metal ’’superalloy” engines. The corrosion properties 
of ceramics make them better candidates to withstand the 
higher fuel impurities expected in the future. 

At present, metal engines require expensive, imported 
metals such as chromium, cobalt, columbium, and nickel. 
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These metals are in limited supply and come from areas of 
potential political instability. Ceramics rely on inexpen- 
sive raw materials which are found in the United States. 

B. DESIGN CONSIDERATIONS 

Ceramics are brittle materials. This does not mean 
that they are necessarily fragile. A brittle material will 
fracture with little or no plastic deformation. A fragile 
material will fracture at low stress with little or no 
plastic deformation. It has been suggested that the word 
"inductile" be used instead of brittle [1] . The inductile 
nature of ceramics is the crux of the design problem with 
this material. Local stress concentrators, such as the tip 
of a crack, will not yield in plastic deformation. If the 
stress concentration is not relieved, material fracture may 
result . 

Another major problem is the variability of data on the 
material properties of ceramics. 

These two problems have a major impact on the design. 
Design procedures that were valid for most engineering 
metals may fail when applied to ceramics. 

C. PROJECT BACKGROUND 

Some of the calculations performed later in this work 
are concerned with the design of a ceramic engine for a 
Defense Advanced Research Projects Agency (DARPA) project. 
The prime contractor, the AiResearch Division of the 
Garrett Corporation, used their T76 model metal engine as 
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the basis for the ceramic design. They indicated that by 
replacing certain metal hot section parts, they could im- 
prove engine output power by 40%, and reduce specific fuel 
consumption by 10%, as shown in Fig. 1. An engine cross- 
section (Fig. 2) shows the ceramic hot section. 

As mentioned earlier, the brittle nature of ceramics 
drives the design. In the case of the Garrett engine, a 
compliant layer was added in the contact region of the 
dovetail as an attempt to reduce the contact stresses. 

Blade tip clearance was increased to reduce the possibility 
of impact of the blade on the shroud. In addition, a prob- 
abilistic approach, based on the Weibull distribution, 
was used as part of a design methodology that was aimed at 
predicting probability of success. These are but a few of 
the examples of how ceramics affected the design. 

D. REASON FOR THESIS 

The complex geometry of many of the ceramic engine's 
parts, and the need to know the stresses to a high degree 
of accuracy, indicated that the finite element method 
should be employed. One analysis code used a finite element 
program based on an isoparametric 12-node brick. This 
allowed quadratic variation in one coordinate direction, 
but only linear variation in the other two coordinate direc- 
tions. This introduced an artificial stiffness into the 
structure. This study was made to determine the effects of 
using a totally quadratic isoparametric brick, vice one 
which is quadratic in only one direction. 
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II. THEORY 



Recall from strength of materials: 

M = EI/p 

where M = moment, E = Young's modulus, I = moment of iner- 
tia, and p = radius of curvature. 

Consider a beam made of an isotropic material, v;ith a 
constant cross-section. If the beam is in pure bending, 
the moment is constant across the span. If Young's modulus 
and moment of inertia are also constant, the radius of cur- 
vature is a constant. A curve which has a, constant radius 
of curvature is an arc of a circle. 

Figure 3(a) shows a 4-node plane element. There can 
only be linear variation on each of the sides. If the ele- 
ment is deformed in pure bending, each of the sides remains 
straight. The top and bottom surfaces do not reflect any 
of the curvature expected. This is a very stiff approxima- 
tion. 

A 6-node plane element, shown in Figs. 3(b) and 3(c), 
yields a better approximation to the bending phenomenon, 
depending on the orientation of the one quadratic side. 
Orientation 1 does not reflect any of the bending curvature. 
Orientation 2, however, does. The finite element method 
would represent the deflected curve, which is defined by 3 
nodes, by passing a parabola through the node points. This 
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can give a very close approximation to the circular arc 
predicted by strength of materials. It should be noted, 
however, that orientation 2 is constrained by the same 
"planes remain plane" restriction applied to slender beam 
theory. Thus, depending on the geometry and loading, a 6- 
node plane element with orientation 2 may or may not accu- 
rately reflect bending for a deep beam. 

A totally quadratic element, shown in Fig. 3(d), is 
not constrained by orientation, geometry or loading. An 8- 
node plane element can give a better representation of bend- 
ing. 

For a given material, a stiff er element, such as the 4- 
or 6-node plane element, will deform less than it should. 
Since strain is based on displacement, and stress is based 
on strain, the stiffen structure will predict a lower stress 
than is actually present. This is a potentially dangerous 
design situation. 

The development presented here is based on a two-dimen- 
sional plane element. Similar reasoning can be applied to 
three-dimensional bricks. A 12-node brick, for example, 
has only one quadratic side. Depending on the loading, 
geometry, and orientation of the quadratic side, a 12-node 
brick could approach the correct solution. It will always, 
however, be stiff er than the actual structure. If it is 
oriented incorrectly with respect to the applied moment, it 
could result in seriously incorrect answers for a coarse 
mesh system. In addition, since two of the coordinate 
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directions are stiffer than the third quadratic coordinate 
direction, an artificial anisotropy is introduced into the 
problem. 

A fully quadratic brick does not have the problems men- 
tioned above. Although the discretization process intro- 
duces some artificial stiffness, it is very slight. A 
fully cubic brick would even be a better representation of 
the actual structure. 
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III. BENDING SPECIMEN PROBLEM 



To evaluate the accuracy of a finite element computer 
program, it is often useful to compare its results to experi- 
mental results for a problem that can be verified experimen- 
tally. A 4-point flexure test has been used very often for 
this purpose [3] . This chapter explains how the flexure 
test can be modeled using 12-node and 20-node bricks. The 
results are also discussed. 

A. FOUR POINT FLEXURE TESTS 

Figure 4 shows a schematic of a 4-point flexure rig. 

It is used to provide test results from a specimen in pure 
bending. Figure (5a) shows a free body diagram for the 
specimen. Figures (5b) and (5c) show the shear and moment 
diagrams. The moment diagram shows that the portion of the 
span between the applied forces on the top surface has no 
shear, and maximum constant moment, i.e., pure bending. 

B. MODELING TECHNIQUES 

A standard 0.125 inch x 0.250 inch x 3.0 inch bend 
specimen was used. The specimen was modeled using a vari- 
ably-noded three-dimension brick from the ADINA [4] and 
SAP IV [5] finite element programs. Maximum use was made 
of the symmetry in the problem. Only 25% of the specimen 
needed to be modeled, as shown in Figs. 6a,b,c. Symmetry 
in the spanwise and long tranverse directions was used. 
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A cross-sectional view of the deformed bar, Fig. 7, shows 
the symmetry in the long transverse direction. The mid-span 
cross-section was given a guided end condition (no displace- 
ment in the x-direction ) . The midplane in the long trans- 
verse direction was also given a guided end condition (no 
displacement in the y-direction) . 

Since the model was very regular in shape, automatic 
nodal point coordinate and connectivity generation were ex- 
tensively used. Both were generated in the spanwise direc- 
tion to take maximum advantage of the automatic generation 
features of the codes employed. 

C. FLEXURE CASE STUDY 

A 12-node brick has one quadratic direction. The other 
two directions are linear. When used to model the bending 
specimen, three distinct orientations of the brick can occur. 
Figure 8 shows each of these orientations with respect to 
the bending moment. An analysis was performed using the 
SAP IV program for each of the three orientations, called 
Case 1, Case 2, and Case 3. The applied forces were modeled 
as distributed loads. 

1. Analysis of Results 

The three-dimensional finite elements used to solve 
the problem would give a solution converging to an exact 
three-dimensional elasticity theory solution. However, such 
a solution in closed form does not exist. Therefore, the 
results were compared to the usual strength of materials 
beam theory. Such a theory is not perfect and circumspection 
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had to be used in such a comparison. For instance, beam 
theory does not predict any anticlastic curvature [18], but 
the finite element solution shows some. However, due to 
the boundary conditions and symmetry, beam theory should 
apply to the x-z plane of symmetry with a very good accura- 
cy. The displacements of the neutral axis of this plane 
from X = 0 to X = L were plotted versus the span length. 
Displacements were normalized by dividing them by the mid- 
span (x = 0) displacement predicted by beam theory 
This facilitates comparisons in terms of the percentage of 
the maximum deflection of beam theory. 

Two loadings and materials were used. The first 
specimen was made of steel, and had a 300 lb load. The 
second specimen was hot pressed silicon nitride (HPSN), 
and had a 180 lb load. The mechanical properties of HPSN 
are given in Ref. [6]. 

The deflection curves were calculated using [7] : 

2 
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- (2L^ + 2aL - a^) 

w(L-a) = reaction end moment 
applied force 
Young’s modulus 
moment of inertia 



Figure 9 shows that L is half the distance between supports 
on the lower side of the specimen, and w is the total load 
applied, i.e., 300 lb or 180 lb. The symbol < > denotes a 
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discontinuous function, i.e., for x<a the quantity <x-a> = 0, 
for x>0 <x-a> = x-a. 

The results of these analyses are shown in Figs. 10 
and 11. The Case 1, 2 and 3 curves correspond to three 
orientations of the 12-node brick. Case 1, the orientation 
with the quadratic side along the axis of bending, is the 
closest of the three to the values predicted by beam theory. 
The maximum deflection of the ceramic specimen is approxi- 
mately 96% of the deflection predicted by beam theory. This 
model is slightly stiffer than beam theory. Cases 2 and 3, 
which have the quadratic side perpendicular to the axis of 
bending, are significantly stiffer models than Case 1. The 
maximum deflection predicted by these cases is only 70% of 
the deflection predicted by beam theory. 

In the displacement-based finite element method, 
stress and displacement are related by [8]: 



{e} = [B] {u} 

and {t} = [C] {e} = [C] [B] {u} 



where 



{t} = a vector composed of the six independent 
components of the stress tensor. 

{e} = a vector composed of the six independent 
components of the small strain tensor. 

IB] = strain-displacement transformation matrix, 
{u} = displacement matrix. 

IC] = matrix of material properties. 
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This formulation clearly shows that a stiffer struc- 
ture (i.e., less displacement and displacement gradients) 
will produce less stress. 

Case 4, shows results which are slightly more flex- 
ible than beam theory. Since beam theory is inherently 
stiff due to the approximations made during its development 
(i.e., negligible shear stress), the Case 4 model is closer 
to predicting the actual stress distribution than either 
Cases 1, 2, 3 or beam theory. It should also be remembered 
that the actual structure modeled will bend in all three 
directions. Clearly, the only way to fully represent such 
a deformation is to use a fully quadratic element. Any 
element that is less than fully quadratic will be too stiff. 

D. TEMPERATURE GRADIENT CASE STUDIES 

The temperature gradient case studies used the same 
models as the flexure case study. Cases 1, 2, 3 refer to 
the same orientation of the 12-node brick as in the flexure 
study. Case 4 is the 20-node brick model. 

Hot pressed silicon nitride was used as the specimen 
material. The coefficient of thermal expansion is given in 
Ref. [6]. A 10 0*^F temperature gradient was imposed across 
the short transverse (0.125 inch). The stress-free refer- 
ence temperature was 1500°F. A variably-noded (12-21) 
three-dimensional brick from the ADINA finite element pro- 
gram was used to model the problem. Nodal temperatures were 
input via a temperature tape. This procedure is described 
in Appendix A. 
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As in the flexure case study, the normalized deflections 
of the y-z plane of symmetry were plotted versus distance 
along the span. Beam theory was used as the basis for com- 
parison. The deflection predicted by beam theory is [7] : 



6 




T^) 



where: 6^ = -aL^(T 2 - T^)/2t 

a = coefficient of thermal expansion ( inch/inch/ °F) 
t = depth of beam 

(T 2 “ = temperature difference through depth of beam 

L = half the distance between the lower reactions 

Two different case studies evolved from the temperature 
problem. Reference [9] states that the Bernoulli-Euler as- 
sumptions can be used in practical analyses of beams under 
thermal loadings. The assumptions state, in part, that the 
effects of lateral contraction can be ignored, i.e., Poisson's 
ratio (v) may be taken as equal to zero. In the case studies 
that follow, Poisson's ratio was used in the first study, and 
set equal to zero in the second. 

1 . Results of the Case Study using Poisson's ration (v) 
The normalized deflection curves for this case study 
are shown in Fig. 12. The displacements from x=0 to x=L are 
plotted. Cases 2 and 3 are stiffer than beam theory. The 
curves show that the maximum deflection calculated for these 
two cases is 75% of the maximum deflection predicted by beam 
theory. 
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The results of Case 4, which used the fully quadratic 
2.0-node brick, show excellent agreement with simple theory. 

Figure 12 shows that the results of Case 1 (the 12- 
node brick with quadratic side along the span) predict a de- 
flection greater than theory. It is puzzling that the re- 
sults show that the element is even more flexible than the 
20-node brick. The temperature algorithm used by ADINA was 
checked extensively. No errors were detected. No satis- 
factory explanation for the 12-node brick behavior has been 
found. It is suspected that the inherent anisotropy of the 
12-node brick may be the cause of the problem. However, 
extensive further examination would be necessary to resolve 
the dilemma. 

Possibly more interesting than the displacement re- 
sults are the stresses calculated. Recall that if a linear 
temperature gradient were imposed on a bar that had no ex- 
ternal constraints, the bar would deform into an arc of a 
circle. Also there would be no stress in the bar according 
to a fundamental theorem of elasticity. In this case study, 
the bottom nodes at x = 0.75 inch were fixed in the z-direc- 
tion. The part of the specimen to the left of this con- 
straint (i.e., toward the free end) should, however, show no 
stress until one reaches the immediate vicinity of the con- 
straint. Each of the cases that use the 12-node brick pre- 
diet stresses, both normal and shear, on the order of 10 

4 

to 10 PSI in the free end. The 20-node brick, however, 
predicts stresses of approximately 5 PSI or less in this 
region. (E=45xl0 psi, v = 0.2, a=1.7xl0 in/in/'^F). 



The above discussion concerning the 12-node brick 
models indicates the possible uncertainty that can result 
from using this element in thermal calculations. 

2 . Results of Case Study with v = 0 

This case study was undertaken in an attempt to ex- 
plain the results of the Case 1 orientation of the 12-node 
brick discussed above. Everything was the same as in the 
previous study, except that the effects of Poisson’s ratio 
were neglected. The results are shown in Fig. 13. The Case 
2 and 3 models have become significantly stiffer. The maxi- 
mum deflection was approximately 67% of the value predicted 
by theory, as compared to approximately 75% when Poisson's 
ratio was used. The results for Case 4 have remained the 
same (to three significant digits). The displacements for 
Case 1 have decreased to such a degree that they were almost 
identical to Case 4. The spurious free end stresses cited 
previously remained. 

In summary, by ignoring the effects of Poisson's 
ratio, three of the models have become significantly stiffer, 
while the fully quadratic model has remained substantially 
the same and in close agreement with the theory. 

It is interesting to note that if the specimen were 
truly unconstrained and had a linear temperature gradient, 
the results of analyses that take into account Poisson's 
ratio should agree exactly with the results of analyses that 
do not 19] . 
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IV. AIRFOIL MODELING 



This chapter discusses the techniques developed to con- 
struct a finite element mesh for the airfoil section of a 
ceramic gas turbine blade. This particular engine airfoil 
lent itself to detailed analysis because of the relative 
ease with which the mesh could be generated. This is due 
to the linear fairing between the tip and root sections of 
the airfoil. Although linear fairing was not optimal from 
an aerodynamic point of view, it was desirable from a manu- 
facturing point of view. Its relative ease of fabrication 
caused it to be used in practice. 

A. DEVELOPMENT OF AIRFOIL PROFILES 

The Garrett Corporation provided Mylar cross-sectional 
profiles of the first stage rotor blade of their experimen- 
tal ceramic turbine [10]. These profiles were for various 
radial distances as measured from the axis of rotation 
along the stacking axis of the blade. In addition to an 
outline of the cross-section and its relation to the stack- 
ing axis, a table of x-y coordinates was given for the pro- 
file. These coordinates described 62 points on the profile. 

The 62 profile points were assembled into x and y ma- 
trices. These matrices were used to develop cross-sectional 
plots on the time sharing terminal and later on the Calcomp 
plotter. Figures 14 and 15 show representative cross- 
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sections of the tip and hub sections. The profile points 
are marked by an *x' and are connected by straight line 
segments . 

Since there was linear fairing between the tip and hub, 
it could be easily verified that the coordinates listed on 
the Mylar drawings were given at the same relative positions 
along the foil profile. In other words, the ith point on 
the R=3.1 inch profile corresponded to the ith point on the 
R=4.0 inch profile, where R is the radial distance along the 
stacking axis from the axis of rotation. 

The root profile was then divided into elements. Care 
was taken to develop elements with convex corners. Also, a 
relatively fine mesh was desirable in the leading and trail- 
ing edge areas of the airfoil to enable accurate modeling 
of the anticipated stress concentrations there. The result- 
ing leading and trailing edge elements proved quite trouble- 
some; this will be discussed more fully below. 

The mesh developed in the airfoil root section dictated 
the mesh in the tip section. This was due to the unique 
correspondence of points in the tip and hub sections. The 
tip section was then connected to the hub section by linear 
nodal coordinate generation. Connectivities were generated 
from the leading to the trailing edge. This scheme took 
maximum advantage of the automatic generation features of 
SAP IV and ADINA. 

B. THE NEGATIVE JACOBIAN DETERMINANT 

As mentioned earlier, the leading and trailing edge ele- 
ments were ill-behaved. The leading and trailing edges of 
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this blade are arcs of a circle. When an element with qua- 
dratic sides is used, these elements are prone to re-entrant 
angles, and resultant negative or zero determinants of the 
Jacobian matrix. In other words, the inverse of the Jacobi- 
an matrix does not exist because there is not a unique cor- 
respondence between the natural and local coordinates of the 
element [8]. This will stop the program because the inverse 
of the Jacobian matrix is used in the strain-displacement 
matrix [B] . Recall: 
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Also: [K] = [B]^[C][B] d(VOL) 

[J] ~^ = inverse of the Jacobian matrix 

[K] = stiffness matrix 

{u} = [H] {Uj^} - displacement field 

[H] = matrix of shape functions 

[C] = matrix of elastic properties 

{e} = strain matrix 

r,s,t = natural coordinates 

As the above equations show, the inverse of the Jacobian 
is an essential part of the calculations for the strain ma- 
trix, and the stiffness matrix. Experience has also shown 

_7 

that when the determinant of the Jacobian is small (10 or 
smaller) the displacements at that node may lose almost all 
their numerical significance. 

C. MESH VERIFICATION THROUGH PSAPl 

Kibler [11] implemented a finite element preprocessor 

! 

for mesh verification for the SAP IV and ADINA programs. 

This enables the user to check visually the coordinates and 
geometry of a mesh by generating graphical displays of the 
completed mesh. Oblique orthographic projections are pro- 
duced using the NPS Calcomp model 765 plotter. Other dis- 
plays are possible. 

The PSAPl program is especially useful because of the 
various mesh presentation options in the program. Among 
them are : 

1. Node numbering 

2. Element numbering 



29 



3. Exploded plots 

4. Full rotation about each axis 

5. Sectioned plots based on geometry, 
node number, or element number 

6. Displacement postprocessing using either 
the deformed structure, or displacement 
vectors at the nodes. (See Chapter VII) 

Other options, as well as a User’s Manual, are fully dis- 
cussed in Ref. [11] . 

Figures 16 through 19 show the various mesh presentations, 
for the linear airfoil sections. PSAPl was used extensively 
in the production of the airfoil mesh. The clear graphical 
presentation assures the engineer that he is using a mesh 
free of geometric or connectivity errors. 

The elements with negative Jacobians were separately 
selected from the mesh and greatly magnified. The re-entrant 
angles were clearly visible on these plots as shown in Fig. 

20. Large scale plots of the hub and tip cross-sections for 
these elements were made by hand. The mid-size nodes were 
moved inward toward the geometrical center of these elements. 

Two airfoil meshes were developed. One was composed of 
20-node bricks, and the other of 12-node bricks. The 12-node 
bricks were oriented such that the quadratic side could be 
used to model the airfoil camber. This orientation was the 
same as the one used by Garrett. Since the 12-node brick 
has only one quadratic side, part of the leading and trailing 
edge curves were eliminated. This is shown in Fig. 21. Aside 
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from this difference, the geometry represented by the two 
meshes was identical. The 12-node airfoil is shown in Fig. 
22 . 

The characteristics of each of these meshes are shown 
in Table I. 
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V. AIRFOIL CASE STUDY 



After thoroughly investigating the ceramic bar problem, 
the next step was to analyze linear airfoil section of the 
ceramic gas turbine blade using the ADINA code 12-node and 
20-node bricks. The bar problem showed that significant 
differences occurred as a result of the brick type and 
orientation used. The purpose of the airfoil study was to 
determine the difference, if any, between the 12-node brick 
airfoil and a 20-node brick airfoil. Each airfoil was 
loaded using the actual design point conditions that could 
be expected in the ceramic gas turbine. Two general types 
of loading were investigated: centrifugal and thermal. 

A. CENTRIFUGAL LOAD 

A centrifugal load preprocessor was developed at NFS, 
and is discussed fully in Ref. [12]. The preprocessor can 
be used to generate consistent nodal loads for a given 
angular speed. The rotor design speed of 41,730 RPM was 
used. 

Figures 23, 24 show the deformed foil superposed on the 
undeformed foil. (Chapter VII discusses using PSAPl to plot 
the deformed structure.) The root plane coincides exactly 
in all cases because the nodes there were fully constrained. 
Two general deformations can be seen: first, there is a 

general elongation of the blade in the radial direction, and 
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second, there is an untwisting of the blade. Figures 25-26 
show the deformed 12-node and 20-node foils superposed on 
each other. It can be seen that the 12-node model is stif- 
fen in representing the untwisting deformation. Tables II 
and III show comparisons of stresses and maximum principal 
stresses, as calculated by the 12-node and 20-node bricks, 
at selected nodes in the root section. At each of the tabu- 
lated nodes, the stresses predicted by the 20-node brick are 
higher than those predicted by the 12-node brick. The 
stress contours for the root section are shown in Figs. 29 
and 30. 



B. THERMAL LOAD 

Reference [13] shows the steady state isotherms. For 
this study, the isotherms were modeled as constant for 
various radial distances along the foil. The cited reference 
shows that this is a good assumption in the root section, 
where the isotherms are flat and closely spaced. In the 
upper portion of the foil, however, the isotherms have a 
definite curvature in the radial direction. Figure 31 shows 
the nodal temperatures used in this study. Preparation of 
the nodal temperature tape is discussed fully in Appendix A. 

Figure 32 shows the deformed foil superimposed on the 
original geometry. A general enlargement of the foil in all 
directions is seen. Table IV shows a comparison of thermal 
stresses at selected points. (It should be noted that these 
are not evaluated at nodal points, but rather at the Gauss 
points. ADINA does not give nodal stresses for thermal load- 
ings . ) 
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C. CONCLUSIONS FROM AIRFOIL STUDY 



This comparative study shows how the choice of a finite 
element can affect the results. The 12-node brick, even 
though oriented in the most advantageous way, is stiffen 
than the fully quadratic three-dimensional element. The 
stiffen element produces lower stresses in the root section. 
It is these root section stresses which transmit the tensile 
forces and moments to the blade platform and eventually to 
the critical dovetail area. 
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VI. DISTRIBUTED LOAD PREPROCESSOR 



The version of ADINA installed at NPS does not contain 
a distributed load feature. Finite element texts, such as 
Ref. [14] , often state a consistent allocation of a uniform 
surface load for a rectangular element. The consistent 
load for a curved surface, such as the airfoil studied in 
this report, is more complex. Since a distributed load is 
used to model the aerodynamic pressure on an airfoil, a 
program to calculate consistent pressure loads for curved 
surfaces was necessary. 

A. DEVELOPMENT OF THEORY 

A consistent load is defined as a set of nodal forces 
which produce the same amount of work as the original force. 
Stated mathematically: 



<v.>{up = / 



AREA 



<f^>{u} d(AREA) 



where: <v^> = vector of nodal consistent loads 



{u^} = nodal displacements 
<f > = vector of surface forces 



{u} =[N£]{u^} = displacement field 
[N^] = shape functions 



Consider the general surface shown in Fig. 33. 




<v . > 
1 




( 1 ) 



1 



35 



n = unit normal to the surface (positive outward), 
p = pressure on dA. 

u = displacement vector = u_i + v^ + wk 

For ease of computation, it is desirable to perform the 
integration in the natural coordinate system. The vector 
ndA is transformed to the natural coordinate system by: 

3R 3R 

ndA = (-^ X -^) da d3 

— C7 OC C7 p 

where: R = xi + yj_ + zk = position vector of point on 

surface as shown in Fig. 33. 
da d3 = differential area in natural coordinate 
system. 

a, 3 are dummy variables standing for the appropri- 
ate natural coordinates which define the 
pressure face (See Table V). 



3R 

3a 



3R 

T3 



This is a vector tangent to a line a where 3 
is constant. The components of this vector 
are elements of the Jacobian matrix. 

is similarly defined. 



The pressure and displacement fields can be discretized 
with the same shape functions : 

p = <N^>{p^} (2) 
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where : 



= shape functions 



{p^} = matrix of nodal pressures 




^i 




matrix of nodal displacements 



Let : 

3R 8R 

^ ~ ^2-i ^ 



then : 



ndA*u = <A^ A 2 ^ 3 ^ 




< V > da dg 
w 



(4) 



^1 " ^2, a “’ 3,3 " ^3, a ^2,B 
^2 " ^3, a ^1,3 ■ ^l,a ^3,3 
^3 ■ ^l,a ^2,3 “ ^2, a ^1,3 



Using Eqs (2 through 5), Eq (1) can now be written 
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<N^>{p^}da d 3 
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<v . > 
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ri 



<Xj^X^X3> 



-1 
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-1 



Ni 

0 

0 



0 



0 

0 
L 

0 N 



N. 

1 



<N^>{p^} da d3 (6) 



B. PROGRAM STRUCTURE 



The consistent pressure load program is structured as 
follows ; 



Subroutine TRANS 

- reads element number, coordinates, and connectivity 

- reads the face to which pressure applied (NFACE) 

- reads pressure at each node of NFACE 

- reads number of Gauss points to be used 
in surface integration (2-6) 

- calls CUBAT • 



Subroutine CUBAT 

- establishes vector of weights and sampling 
points to be used in Gaussian integration 

- calls SHAPE 

Subroutine SHAPE 

- evaluates shape functions at Gauss points 

- evaluates derivative of shape function at 
Gauss point 

- forms Jacobian matrix 

- evaluates determinant of Jacobian matrix to check 
for zero or negative Jacobian determinant 

- calls FACEP 
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Subroutine FACEP 



- establishes matrices in integrand 

- performs matrix multiplication to establish 
integrand 

- performs surface integration 

A listing of the consistent pressure load program is in 
Appendix C. 

C. PROGRAM VERIFICATION 

Several test cases were run to insure the accuracy of 
the pressure load preprocessor. A uniform distributed sur- 
face load was rotated through all faces of a one inch cube, 
and the 2-6 Gauss point integration scheme was checked for 
each face. The results agreed exactly, both in magnitude 
and sign, with the consistent load allocation found in 
Ref. [14] . 

The program was successfully checked to insure that dis- 
tributed surface force could be applied to more than one 
face of the element simultaneously, and that pressures could 
be applied simultaneously to faces that share a common edge. 

The ceramic flexure specimen mesh (Case 4) was also 
used to verify an actual problem. Pressure loads were 
applied to various faces of the specimen, including adjoin- 
ing faces. In each case, the program produced correct re- 
sults . 

As a final test, a distributed load was also placed 
across the span between the supports. A sketch of the 
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problem is shown in Fig. 34. The results were again com- 
pared with simple beam theory. The results were excellent 
as shown in Fig. 35. 



D. USER’S MANUAL 



Figure 36 shows a hexahedral element in the natural co- 
ordinate system. The node numbering convention is also 
shown in Fig. 36. Table V shows the face numbering conven- 
tion used in the pressure load preprocessor. Each face has 
four corner nodes and four midside nodes. The pressure dis- 
tribution acting on an element face is defined by specifying 
the nodal pressure intensities at the corner and mid-side 
nodes. For example, if Face 4 were loaded Pg, Py, Pg, 

^20’ ^15’ ^19’ ^11 specified in that order. 

Appendix B lists the input data and results for the sample 
bar problem, shown in Fig. 36. 

1. Control Card (1015) 



Columns 



1-5 



6-10 



11-15 



16-20 



Variable 

NUMNP 



NUMEL 



NGP 



NCUR 



Meaning 

Total number of node 
points in problem. 

Total number of 
elements in problem. 

Number of Gauss points 
to use in surface 
integration. (2<NGP<6) 
Number of ADINA load 
curve that output loads 
refer to. (usually ,NCUR=1) 



40 



21-25 



26-30 



31-35 



36-40 



41-45 



2 . 



1-5 

6-10 



11-50 



NPLOAD 
IPRINT 1 
IPRINT 2 
IPRINT 3 
IPRINT 4 

Nodal Pressures Card (215, 

NPREL (J) 
NFACE (J) 

PRESS (I,J) 



Number of element 
faces with distributed 
loads . 

1-program will just 
read and print ADINA 
input deck 
1-program will not 
print element-by- 
element solution. 
1-program will not 
punch output deck of 
consistent loads. 
1-program will omit 
printing input data. 

8F5. 0) 

element number 
face to which pressure 
applied. (See Table V) 
nodal pressure 
intensities input here 
corner nodes first, 
then midside. (See 
Table V) Positive pres 
sure is directed in 
same sense as outward 
unit normal. 
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E. USE OF PRESSURE LOAD PREPROCESSOR ON CP/CMS 

The pressure load preprocessor is ideally -suited for use 
on CP/CMS. A compiled version of the program (f iletype=text ) 
should be stored on a private disk or workspace. A data 
file, such as FT04F002, should also be defined. The input 
data deck (Appendix B) can then be read into this file [15]. 

It is suggested that an executive type routine also be 
used: 

S TYPEOUT ERROR 

FILEDEF 05 DSK FILE FT04F002 (PERM) 

FILEDEF 06 PRT (PERM) 

FILEDEF 07 PUN (PERM) 

VSET RDYMSG OFF 
$ PRESSURE 

This defines your input/output files and executes the 
program. Output is directed to the line printer, and card 
punch . 

The pressure program utilizes approximately 3 0 OK BYTES 
of core. The extra core should be requested at LOGIN by: 
L/^nnnnPtt/^ 3 0 OK 

where : 

nnnn = user number 
tt = terminal number 



42 



VII. USING PSAPl AS A DISPLACEMENT POSTPROCESSOR 



PSAPl is capable of postprocessing displacements. The 
results can be presented in two ways: 

1. The deformed structure can be plotted. In this case, 
scaled displacements are added to the original geometry to 
produce the deformed coordinates. The deformed structure is 
then plotted using this new geometry. 

2. The displacements can be represented as scaled vec- 
tors. The scaled displacement vectors are plotted at origi- 
nal geometry nodes . 

The purpose of this chapter is to explain the modifica- 
tion of the ADINA code to store the displacements, and how 
to use the PSAPl postprocessing capabilities. 

A. MODIFICATION OF ADINA 

Reference [16] states that a version of SAP IV had been 
modified to punch a deck of displacement cards from the out- 
put data. When dealing with large quantities of data, how- 
ever, storage of this data on direct access storage devices 
is preferable. Therefore, the ADINA code was modified to 
store the output displacements on device 58, in addition to 
the usual printed output. 

The flag for saving the displacements is the integer 
number 1 in column 80 of the ADINA master control card [4] . 
The flag variable name is ITP 58, and is passed in common 
/EAST/. 
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The job control language (JCL) for ADINA must also be 

modified to define device 58: 

//FT58F001 DD DISP=OLD ,UNIT= 3 3 30 , VOL=SER=DISKO , 

// DSN=S nnnn . pppp 

where : 

nnnn is the user number 

pppp is the name assigned by the user to identify 
the data set. 

Since an unformatted write statement is used to store 
the data on device 58, a variable blocksize must be used 
when allocating disk space. Reference [17] explains this 
procedure . 

B. MODIFICATIONS OF PSAPl 

There are presently two subroutines in PSAPl that can 
be used to input displacement data to PSAPl. 

1. Subroutine DATA 9 was added to PSAPl by Losh [16] . 
This subroutine reads displacement data from punched cards. 
Reference [16] discusses its use, and should be consulted. 

2. Subroutine DATA 5 reads displacement data directly 
from device 58. A JCL card must be added to the PSAPl JCL 
to define device 58: 

//GO.FT58001 DD DISP=0LD ,UNIT= 33 3 0 , V0L=SER=DISK0 , 

// DSN=Snnnn . pppp 

C. USER'S MANUAL FOR DISPLACEMENT POSTPROCESSING 

The PSAPl User's Manual [11] applies with the following 
modifications : 

a) Namelist Option 
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Variable, Value 



Description 



NUDISP = 1 


X direction displacements input 


NVDISP = 1 


y direction displacements input 


NWDISP = 1 


z direction displacements input 


KDATA = 5 


subroutine DATA 5 used to read 
displacement data from device 58 


NVALUS = nnn 


nnn is the integer number of dis 
placements set to be plotted. 
This is equal to the number of 
node points . 



b) Namelist Piet 



KDISP = n 


n=l: plot of deformed structure 

n=3: displacements represented 

by vectors at no^es. 


IDMAG = n 


n=l: direct magnification of 

displacement data by DMAGS 
n=2: scaling of displacement 

data to a maximum value 
of DMAGS. 


DIIAGS = r 


r is a floating point number 
which indicates the magnifica- 
tion of the displacements. (0.1 
has been found appropriate.) 


Figures 37 and 38 


are examples of displacement postpro- 



cessing. The structure is the bending specimen loaded in 
four-point flexure. 
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Figure 37 illustrates the deformed structure option of 
PSAPl. Figure 38 shows the displacements of element 15 
plotted as scaled vectors at the node points. 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

1. The 12-node brick is a stiffer model of a structure 
than the fully quadratic 20-node brick. 

2. The accuracy of the results obtained using a 12-node 
brick depends on geometry, loading, and the orientation of 
the element. 

3 ; The flexure problem showed that two orientations of 
the 12-node brick are stiffer than the third. Resulting dis- 
placements and stresses differed significantly, and in all 
cases were lower than the 20-node brick and theory. 

4. The airfoil study showed that the maximum principal 
stresses predicted in the root section could differ by as 
much as 14.9% depending on the basic element used. The 12- 
node brick under-predicted the 20-node results. The use of 
the 12-node brick could therefore produce an unsafe design. 

B. RECOMMENDATIONS FOR FUTURE VJORK 

As is true with many projects, this thesis raised many 
new questions. The following areas of future work are 
suggested : 

1. The ADINA code should be reviewed with an eye toward 
streamlining it. A basic code using only two- and three- 
dimensional elements, and the linear isotropic material 
model would be a useful tool. A reduced version of ADINA 



covld then be installed on the time-sharing system, or 
other stand-alone system. 

2. The Field program, developed at NFS, should be ex- 
amined and expanded if necessary. This program would serve 
as an excellent thermal preprocessor. It should be modi- 
fied to produce a nodal temperature tape that could be 
used in ADINA, or other code capable of using temperature 
inputs . 

3. Aerodynamic pressure loads should be examined. A 
program that could calculate the nodal pressure intensi- 
ties for an element would provide the necessary input in- 
formation for the distributed load preprocessor. 

4. Centrifugal, pressure and thermal loadings should 
be superposed. 
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SPECIFIC FUEL CONSUMPTION, LB/HP-HR 




OUTPUT SHAFT HORSEPOWER 



Figure 1 



Engine Design Parameters 
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Figure 2, Ceramic Hot-Section Components 




a. 4 Node Plane Element 





c. 6 Node Plane Element: Orientation 2. 




Figure 3. Two-Dimensional Plane Element 
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SCHEMATIC SHOWING SELF-ALIGNING CAPABILITY 



Figure 4. Schematic of 4-Point 
Flexure Rig 



52 



F F 



1.125" 


0.75" 


1 1.125" 




f \ 


f 


/ 

0.75” 


^ / 

- 1 , - 


\ 

0.75" 


A* 







Figure 5a. Bend Specimen Free-Body Diagram 



F 



F 



Figure 5b. 



Bend Specimen Shear Diagram 




Figure 5c. Bend Specimen Moment Diagram 
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Figure 6a. Four-point Bend Specimen Model with Midspan 
Boundary Condition. 
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Figure 6b. x-z View of Half-Span with Guided 
End Condition 
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Figure 6c. y-z View of Specimen with Guided 
Midplane Boundary Condition 
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PLANE OF 
SYMMETRY 



Figure 7. Cross-Sectional View of Deformed 
Specimen Showing x-z Plane of 
Symmetry 
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Figure 8a. 12-Node Brick Orientation for Case-1 




Figure 8b. 12-Node Brick Orientation for Case-2 




Figure 8c. 12-Node Brick Orientation for Case-3 
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Figure 9a. Free Body Diagram for Flexure Case 
Study 




Figure 9b. Free Body Diagram for Temperature Gradient 
Case. Study 
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Figu:'>e 10. Normalized Displacement vs. Span for Steel 
Bend Specimen 
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Figure 11. Normalized Displacement vs. Span for Ceramic 
Bend Specimen 
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Figure 12. Normalized Displacement vs. Span for Temperature 
Gradient Analyses (v=0.2 ) 
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Figure 13. Normalized Displacement vs. Span for 
Temperature Gradient Analyses (v=0) 
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Figure 14. Airfoil Tip Section Profile 
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Figure 15. Airfoil Hub Section Profile 
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Figure 16. 20-Node Airfoil Mesh View 




66 



Figure 17. 20-Node Airfoil Mesh View 
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Elements Numbered 




Figure 20. Trailing Edge Element Showing 
Re-entrant Angle 
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Figure 21. 12-Node Airfoil Section Showing 

Truncation of Leading Edge Element 
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Figure 22, 12-Node Airfoil Mesh View 
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of Nodal Points Used for Comparison in 
Tables II and III 
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igure 31, Nodal Temperature Distribution 
used for Thermal Analysis 
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Figure 32. Thermally Deformed 20-Node Airfoil 
Superposed on Original Geometry 
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Figure 34. Pressure Load Verification Problem 
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Figure 35. Normalized Displacement vs. Span for Bar 
Problem with Pressure Load 
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Figure 36. Hexahedral Element in Natural 
Coordinates 
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Figure' 38. Postprocessing Using Displacement Vectors 
at the Nodes (Flexure Problem) 



TABLE I: MESH CHARACTERISTICS 



20-NODE BRICKS 

Element type 20-node brick 

Number of nodes 452 

Number of elements 5 2 

Degrees of freedom 1152 

Mean half -bandwidth 9 3 

Maximum half-bandwidth 12 9 

Number of stiffness matrix elements 107064 



12-NODE BRICKS 



Element type 12-node brick 

Number of nodes 27 0 

Number of elements 52 

Degrees of freedom 648 

Mean half-bandwidth 4 5 

Maximum half-bandwidth 6 6 

Number of stiffness matrix elements 29106 
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TABLE II. COMPARISON OF MAXIMUM STRESSES IN AIRFOIL ROOT 
SECTION (Cf. Fig. 27) 



LOCATION 


FOIL 


a 

XX 

(PSD 


cr 

yy 

(PSD 


(PSD 


A 


12 


2553 


2553 


10213 




20 


2658 


2658 


10634 


B 


12 


2897 


2897 


11588 




20 


3211 


3211 


12843 


C 


12 


3082 


3082 


12329 




20 


3560 


3560 


14244 


D 


12 


3130 


3130 


12520 




20 


3656 


3656 


14626 


E 


12 


3074 


3074 


12300 




20 


3625 


3625 


15000 


F 


12 


2946 


2946 


11784 




20 


3447 


3447 


137 8 8 


6 


12 


2781 


2781 


11125 




20 


3231 


3231 


12924 


H 


12 


2608 


2608 


10432 




20 


3000 


3000 


12001 
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TABLE III. COMPARISON OF MAXIMUM PRINCIPAL STRESSES IN 
AIRFOIL ROOT SECTION (cf. Fig. 27) 



LOCATION 


12-NODE AIRFOIL 
(PSD 


*^1 

20-NODE AIRFOIL 
(PSD 


% DIFFERENCE 


A 


10441 


10865 


3.9 


B 


11790 


12992 


9.25 


C 


12519 


14379 


12.9 


D 


12706 


14759 


13.9 


E 


12486 


14676 


14.9 


F 


11969 


13998 


14.5 


G 


11308 


13192 


14,3 


H 


10608 


12293 


13.7 



a^(20) - a (12) 

% Difference E 

a^( 20 ) 
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TABLE IV. COMPARISON OF SELECTED THERMAL STRESSES IN 
AIRFOIL ROOT SECTION (cf. Fig. 28) 



LOCATION 


^zz 

12 NODE 
(PSD 


a 

zz 

20 NODE 
(PSD 


% DIFF 


A 


-11590 


-11407 


-1.6 


B 


-7744 


-8664 


10.6 


C 


-4387 


-5759 


23.8 


D 


-1813 


-3487 


48.0 


E 


-364 


-2337 


84.0 


F 


169 


-1800 


109.4 


G 


-668 


-2431 


72.5 


H 


-3560 


-4998 


28.8 


I 


-10543 


-12554 


16.01 



% Difference 



a(20) - a(12) 
a(20) 



91 



TABLE V. FACE NUMBERING CONVENTION FOR PRESSURE LOADS 
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APPENDIX A: PREPARATION OF TEMPERATURE TAPE 



If temperature input is used with the ADINA program, 
nodal point temperatures are stored in sequence on device 
number 56 [4], At NPS, ITEL 3330 disks have been used, 
although any storage device may be designated. NSTE + 1 
records are read, where NSTE is the number of solution 
steps. (In linear static analysis, NSTE usually equals 1.) 

The read statement in the program is: 

READ (56) TIME, (TEMPV(I), 1=1, NUMNP). 

On each record, the problem time of the nodal tempera- 
tures is read, followed by an implicit "do-loop" from 1 to 
the total number of node points (NUMNP) to read each nodal 
point temperature. 

A program to establish the temperature tape is listed 
at the end of this Appendix. It will be discussed in three 
parts; (A) job control language, (B) the RDADIN subroutine 
and (C) the binary write. 

A. TEMPERATURE TAPE JOB CONTROL LANGUAGE (JCL) 

The program listing in Appendix B starts with a two step 
procedure [17] . The first is a purge operation that returns 
the previously reserved disk space to the operating system, 
the second step reserves disk space for the data set which 
is named F0099.TAPE. One cylinder of 3330 disk space (248K 
bytes) is reserved on one of the NPS computer center resident 
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disks. The data set has a variable block size, and an ex- 
piration date of 360 days after establishment. 

The JCL employed combines ease of continued use of the 
disk space, with minimum chances of data set expiration. 

B. THE RDADIN SUBROUTINE 

The RDADIN subroutine was originally programmed as part 
of the PSAPl preprocessor. The subroutine will read an 
ADINA input deck from title card through element connectivity. 
It generates node numbers, nodal coordinates, element numbers 
and element connectivities for the entire mesh. Each of the 
above is stored in a linear array. The programmer is then 
able to utilize any or all of these arrays to input the 
nodal point temperature. 

For example, in the linear airfoil mesh there were iso- 
thermal planes. The matrix of z-coordinates (radial distance 
from axis of rotation) was used. Each z-coordinate was 
checked by a series of *if statements*, and nodal tempera- 
tures were assigned. 

C. THE BINARY WRITE 

After the matrix of nodal temperatures was established, 
the next step was to write it sequentially onto device 56. 
Array RCDl contained the problem time for the temperatures, 
and temperatures for each node; i.e., it was a linear array 
of dimension NUMNP + 1. This array could then be written 
onto device 56 by the statement: 

WRITE (56)RCD1. 
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Device 56 must then be defined in the GO step by: 

//GO.FT56F001 DD UNIT= 3 3 3 0 , VOL=SER-DISK01 , 

///vDSN=F0 09 9.TAPE,DISP = OLD 

Each time the WRITE statement is executed one record is 
established. The sequence must be repeated until NSTE + 1 
records are established, where NSTE is the number of solu- 
tion steps. 

In the linear airfoil problem, there was only one solu- 
tion time step. Therefore, only two records on device 56 
were needed. One record was at time=0, and the next record 
was one time step increment later. 

It is recommended that a hard copy of the RCDl matrix 
be printed to ascertain the accuracy of the temperature 
tape . 
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APPENDIX B. NODAL TEMPERATURE TAPE PROGRAM 
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APPENDIX C. LISTING OF INPUT DECK AND RESULTS FOR 
SAMPLE PRESSURE PROBLEM 
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NEC ANC NCARO FOR IC IN GEOMl = 642 
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THE TOTAL NUMBER OF LOADS IS 



APPENDIX D. LISTING OF DISTRIBUTED LOAD PREPROCESSOR 
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